package scu.maqiang.fes;

import static scu.maqiang.numeric.Constants.GP2;
import static scu.maqiang.numeric.Constants.GW2;
import static scu.maqiang.numeric.MVO.*;
import static scu.maqiang.numeric.MVO.matmul;

import java.util.Arrays;

import scu.maqiang.mesh.Mesh;
import scu.maqiang.mesh.ScalarFunc;
import scu.maqiang.numeric.CG;
import scu.maqiang.numeric.Constants;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.Pair;
import scu.maqiang.numeric.ParamCheck;
import scu.maqiang.numeric.SRMatrix;

public class FES3H81 extends FES {

	public FES3H81(Mesh mesh) {
		super(mesh, 1);
		eleStiff = FES3H81::elementStiff;
		eleMass = FES3H81::elementMass;
		eleSource = FES3H81::elementSource;
		eleFlux = FES3T41::elementHeatFlux;
	}

	public double getValue(double[] u, double[] xyz) {
		ParamCheck.checkEqual(u.length, this.ndof);
		int[] eleIdx = new int[this.nPerEle];
		double[] eleVal = new double[this.nPerEle];
	    double[] uu = new double[8];
	    double value = 0;
	    double[][] coord = new double[8][3];
	    for(int i = 0, n = mesh.getNt(); i < n; i++) {
	    	int[] ele = mesh.getElement(i);
	    	mesh.getCoordPerElement(ele, coord);
	    	if (CG.pointInH8(coord, xyz)) {
	    		CG.H8ShapeFunction(coord, xyz, eleVal);
	    		for(int j = 0; j < 8; j++) {
	    			uu[j] = u[ele[j]];
	    			eleIdx[j]= ele[j];
	    			value += uu[j] * eleVal[j];
	    		}
	    		return value;
	    	}	    	
	    }
	    throw new IllegalArgumentException("The Point ( " + Arrays.toString(xyz) + " ) is not in the domain!");
	}
	
	public double getValue(double[] u, double[] xyz, int[] eleIdx, double[] eleVal) {
		ParamCheck.checkEqual(u.length, this.ndof);
		ParamCheck.checkEqual(eleIdx.length, this.nPerEle);
		ParamCheck.checkEqual(eleVal.length, this.nPerEle);
	    double[] uu = new double[8];
	    double value = 0;
	    double[][] coord = new double[8][3];
	    for(int i = 0, n = mesh.getNt(); i < n; i++) {
	    	int[] ele = mesh.getElement(i);
	    	mesh.getCoordPerElement(ele, coord);
	    	if (CG.pointInH8(coord, xyz)) {
	    		CG.H8ShapeFunction(coord, xyz, eleVal);
	    		for(int j = 0; j < 8; j++) {
	    			uu[j] = u[ele[j]];
	    			eleIdx[j]= ele[j];
	    			value += uu[j] * eleVal[j];
	    		}
	    		return value;
	    	}	    	
	    }
	    throw new IllegalArgumentException("The Point ( " + Arrays.toString(xyz) + " ) is not in the domain!");
	}
	
	public Pair<Double, Integer> getValue(double[] u, double[] xyz, int previousEleIdx, double[][] coef, int[] eleIdx, double[] eleVal) {
		ParamCheck.checkEqual(u.length, this.ndof);
		ParamCheck.checkEqual(eleIdx.length, this.nPerEle);
		ParamCheck.checkEqual(eleVal.length, this.nPerEle);
		double value = 0;
		double x = xyz[0], y = xyz[1], z = xyz[2];
		double[] v = {1.0, x, y, z, x * y, y * z, z * x, x * y * z};
		double[][] coord = new double[8][3];
		int[] ele = mesh.getElement(previousEleIdx);
		mesh.getCoordPerElement(ele, coord);
		if (CG.pointInH8(coord, xyz)) {
			MVO.matmul(coef, v, eleVal);
			//CG.H8ShapeFunction(coord, xyz, eleVal);
			for (int j = 0; j < 8; j++) {
				//eleIdx[j] = ele.idx[j];
				value += u[ele[j]] * eleVal[j];
			}
			return new Pair<Double, Integer>(value, previousEleIdx);
		} else {
			for (int i = 0, n = mesh.getNt(); i < n; i++) {
				ele = mesh.getElement(i);
				mesh.getCoordPerElement(ele, coord);
				if (CG.pointInH8(coord, xyz)) {
					//CG.H8ShapeFunction(coord, xyz, eleVal);
					CG.H8ShapeFunctionCoef(coord, coef);
					MVO.matmul(coef, v, eleVal);
					for (int j = 0; j < 8; j++) {
						eleIdx[j] = ele[j];
						value += u[ele[j]] * eleVal[j];
					}
					return new Pair<Double, Integer>(value, i);
				}
			}
			throw new IllegalArgumentException("The Point ( " + Arrays.toString(xyz) + " ) is not in the domain!");
		}

	}
	
	public  void getValue(double[] u, double[] xyz, Pair<Double, Integer> data, double[][] coef, int[] eleIdx, double[] eleVal) {
		ParamCheck.checkEqual(u.length, this.ndof);
		ParamCheck.checkEqual(eleIdx.length, this.nPerEle);
		ParamCheck.checkEqual(eleVal.length, this.nPerEle);
		double value = 0;
		double x = xyz[0], y = xyz[1], z = xyz[2];
		double[] v = {1.0, x, y, z, x * y, y * z, z * x, x * y * z};
		double[][] coord = new double[8][3];
		int[] ele = mesh.getElement(data.getSecond());
		mesh.getCoordPerElement(ele, coord);
		if (CG.pointInH8(coord, xyz)) {
			MVO.matmul(coef, v, eleVal);
			//CG.H8ShapeFunction(coord, xyz, eleVal);
			for (int j = 0; j < 8; j++) {
				//eleIdx[j] = ele.idx[j];
				value += u[ele[j]] * eleVal[j];
			}
			data.setFirst(value);
		} else {
			for (int i = 0, n = mesh.getNt(); i < n; i++) {
				ele = mesh.getElement(i);
				mesh.getCoordPerElement(ele, coord);
				if (CG.pointInH8(coord, xyz)) {
					//CG.H8ShapeFunction(coord, xyz, eleVal);
					CG.H8ShapeFunctionCoef(coord, coef);
					MVO.matmul(coef, v, eleVal);
					for (int j = 0; j < 8; j++) {
						eleIdx[j] = ele[j];
						value += u[ele[j]] * eleVal[j];
					}
					data.setFirst(value);
					data.setSecond(i);
					return;
				}
			}
			throw new IllegalArgumentException("The Point ( " + Arrays.toString(xyz) + " ) is not in the domain!");
		}

	}
	
	public  void getValue_Cube(double[] u, double[] xyz, Pair<Double, Integer> data, double[][] coef, int[] eleIdx, double[] eleVal) {
		ParamCheck.checkEqual(u.length, this.ndof);
		ParamCheck.checkEqual(eleIdx.length, this.nPerEle);
		ParamCheck.checkEqual(eleVal.length, this.nPerEle);
		double value = 0;
		double x = xyz[0], y = xyz[1], z = xyz[2];
		double[] v = {1.0, x, y, z, x * y, y * z, z * x, x * y * z};
		double[][] coord = new double[8][3];
		int[] ele = mesh.getElement(data.getSecond());
		mesh.getCoordPerElement(ele, coord);
		if (CG.pointInCube(coord, xyz)) {
			MVO.matmul(coef, v, eleVal);
			//CG.H8ShapeFunction(coord, xyz, eleVal);
			for (int j = 0; j < 8; j++) {
				//eleIdx[j] = ele.idx[j];
				value += u[ele[j]] * eleVal[j];
			}
			data.setFirst(value);
		} else {
			for (int i = 0, n = mesh.getNt(); i < n; i++) {
				ele = mesh.getElement(i);
				mesh.getCoordPerElement(ele, coord);
				if (CG.pointInCube(coord, xyz)) {
					//CG.H8ShapeFunction(coord, xyz, eleVal);
					CG.H8ShapeFunctionCoef(coord, coef);
					MVO.matmul(coef, v, eleVal);
					for (int j = 0; j < 8; j++) {
						eleIdx[j] = ele[j];
						value += u[ele[j]] * eleVal[j];
					}
					data.setFirst(value);
					data.setSecond(i);
					return;
				}
			}
			throw new IllegalArgumentException("The Point ( " + Arrays.toString(xyz) + " ) is not in the domain!");
		}
		//System.out.println(MVO.toString(eleMatrix));
	}

	/**
	 * 组装热刚度矩阵, 材料系数为常数
	 *
	 * @param matrixCoef
	 *            常系数
	 * @param tp
	 *            组装刚度矩阵类型
	 * @param A
	 *            组装得到的稀疏矩阵
	 */
	public void assembleHeatStiff(double[][][] matrixCoef, BVPType tp, SRMatrix A, int... domainLabel) {
		assembleGlobalMatrix(matrixCoef, FES3H81::elementStiff, tp, A, domainLabel);
	}

	/**
	 * ���㵥Ԫ�նȾ���
	 * @param ele
	 * @param xyz
	 * @param kxx
	 * @param tp
	 * @param eleMatrix
	 */
	public static void elementStiff(int[] ele, double[][] xyz, double[][] kxx, BVPType tp, double[][] eleMatrix) {
		double[] x = MVO.col(xyz, 0);
		double[] y = MVO.col(xyz, 1);
		double[] z = MVO.col(xyz, 2);
		double[] dx = new double[8];
		double[] dy = new double[8];
		double[] dz = new double[8];
		double[] v = new double[8];
		double xi, et, zt, detJ, wxarea, c;
		MVO.fill(eleMatrix, 0.0);
		double[][] B = new double[3][8];
		double[][] KK = new double[8][8];
		switch (tp) {
		case COMMON:
			for (int i = 0; i < 2; i++) {
				xi = GP2[i];
				for (int j = 0; j < 2; j++) {
					et = GP2[j];
					for (int k = 0; k < 2; k++) {
						zt = GP2[k];
						detJ = CG.ISOMap3D(CG::H8ShapeFunction, x, y, z, xi, et, zt, v, dx, dy, dz);
						c = MVO.dot_product(kxx[0], v);
						wxarea = GW2[i] * GW2[j] * GW2[k] * detJ * c;
						for (int si = 0; si < 8; si++) {
							for (int sj = 0; sj < 8; sj++) {
								eleMatrix[si][sj] += (dx[si] * dx[sj] + dy[si] * dy[sj] + dz[si] * dz[sj]) * wxarea;
							}
						}
					}
				}
			}
			break;
		case CONSTITUTE_MATRIX_COMMON:
			for (int i = 0; i < 2; i++) {
				xi = GP2[i];
				for (int j = 0; j < 2; j++) {
					et = GP2[j];
					for (int k = 0; k < 2; k++) {
						zt = GP2[k];
						detJ = CG.ISOMap3D(CG::H8ShapeFunction, x, y, z, xi, et, zt, v, dx, dy, dz);
						wxarea = GW2[i] * GW2[j] * GW2[k] * detJ;
						B[0] = dx;
						B[1] = dy;
						B[2] = dz;
						mul(matmul(transpose(B), matmul(kxx, B)), wxarea, KK);
						for (int si = 0; si < 8; si++) {
							for (int sj = 0; sj < 8; sj++) {
								eleMatrix[si][sj] += KK[si][sj];
							}
						}
					}
				}
			}
			break;
		default:
			throw new IllegalArgumentException("Wrong BVPType!");
		}
//		System.out.println(MVO.toString(eleMatrix));
//		System.exit(1);
	}

	public static void elementMass(int[] ele, double[][] xyz, double[][] coef, BVPType tp, double[][] eleMatrix) {
		double[] x = MVO.col(xyz, 0);
		double[] y = MVO.col(xyz, 1);
		double[] z = MVO.col(xyz, 2);
		double[] dx = new double[8];
		double[] dy = new double[8];
		double[] dz = new double[8];
		double[] v = new double[8];
		double xi, et, zt, detJ, wxarea, c;
		MVO.fill(eleMatrix, 0.0);
		switch (tp) {
		case CONSIST_MASS:
			for (int i = 0; i < 2; i++) {
				xi = GP2[i];
				for (int j = 0; j < 2; j++) {
					et = GP2[j];
					for (int k = 0; k < 2; k++) {
						zt = GP2[k];
						detJ = CG.ISOMap3D(CG::H8ShapeFunction, x, y, z, xi, et, zt, v, dx, dy, dz);
						c = MVO.dot_product(coef[0], v);
						wxarea = GW2[i] * GW2[j] * GW2[k] * detJ * c;
						for (int si = 0; si < 8; si++) {
							for (int sj = 0; sj < 8; sj++) {
								eleMatrix[si][sj] += v[si] * v[sj] * wxarea;
							}
						}
					}

				}
			}
			break;
		case LUMP_MASS:
			for (int i = 0; i < 2; i++) {
				xi = GP2[i];
				for (int j = 0; j < 2; j++) {
					et = GP2[j];
					for (int k = 0; k < 2; k++) {
						zt = GP2[k];
						detJ = CG.ISOMap3D(CG::H8ShapeFunction, x, y, z, xi, et, zt, v, dx, dy, dz);
						c = MVO.dot_product(coef[0], v);
						wxarea = GW2[i] * GW2[j] * GW2[k] * detJ * c;
						for (int si = 0; si < 8; si++) {
							eleMatrix[si][si] += v[si] * wxarea;
						}
					}

				}
			}
			break;
		default:
			break;
		}
	}
	
	public void assembleHeatThetaScheme(ScalarFunc[] coefFunc, double[] param, BVPType[] tp, SRMatrix LMatrix, SRMatrix RMatrix) {
		double[][] kxx = new double[1][nPerEle];
		double[][] rho = new double[1][nPerEle];
		double[][] eleMatrixStiff = new double[dofPerElement][dofPerElement];
		double[][] eleMatrixMass = new double[dofPerElement][dofPerElement];
		double[][] eleLMatrix = new double[dofPerElement][dofPerElement];
		double[][] eleRMatrix = new double[dofPerElement][dofPerElement];
		double dt = param[0];
		double theta = param[1];
		for (int i = 0; i < nE; i++) {
			int[] ele = mesh.getElement(i);
			mesh.getCoordPerElement(ele, coordRegion);
			int label = mesh.getElementLabel(i);
			for (int j = 0, nj = nPerEle; j < nj; j++) {
				kxx[0][j] = coefFunc[0].action(coordRegion[j], label, param);
				rho[0][j] = coefFunc[1].action(coordRegion[j], label, param);
			}
			elementStiff(ele, coordRegion, kxx, tp[0], eleMatrixStiff);
			elementMass(ele, coordRegion, rho, tp[1], eleMatrixMass);
			for(int ii = 0; ii < dofPerElement; ii++ ) {
				for(int jj = 0; jj < dofPerElement; jj++) {
					eleLMatrix[ii][jj] = eleMatrixMass[ii][jj] + theta * dt * eleMatrixStiff[ii][jj];
					eleRMatrix[ii][jj] = eleMatrixMass[ii][jj] - (1 - theta) * dt * eleMatrixStiff[ii][jj];
				}
			}
			LMatrix.assemble(ele, eleLMatrix);
			RMatrix.assemble(ele, eleRMatrix);
		}
	}

	public static void elementSource(int[] ele, double[][] xyz, double[][] f, BVPType tp, double[] eleF) {
		double[] x = MVO.col(xyz, 0);
		double[] y = MVO.col(xyz, 1);
		double[] z = MVO.col(xyz, 2);
		double[] dx = new double[8];
		double[] dy = new double[8];
		double[] dz = new double[8];
		double[] v = new double[8];
		double xi, et, zt, detJ, wxarea, c;
		MVO.setValue(eleF, 0.0);
		switch (tp) {
		case COMMON:
			for (int i = 0; i < 2; i++) {
				xi = GP2[i];
				for (int j = 0; j < 2; j++) {
					et = GP2[j];
					for (int k = 0; k < 2; k++) {
						zt = GP2[k];
						detJ = CG.ISOMap3D(CG::H8ShapeFunction, x, y, z, xi, et, zt, v, dx, dy, dz);
						c = MVO.dot_product(f[0], v);
						wxarea = GW2[i] * GW2[j] * GW2[k] * detJ * c;
						for (int si = 0; si < 8; si++) {
							eleF[si] += v[si] * wxarea;
						}
					}
				}
			}
			break;
		default:
			break;
		}
//		System.out.println(Arrays.toString(eleF));
	}
	
	public void assembleHeatConvectSurface(ScalarFunc[] coefFunc, double[] param, BVPType tp, SRMatrix A, int... bdLabel) {
		assembleBoundaryMatrix(coefFunc, param, FES3H81::elementConvectSurface, tp, A, bdLabel);
	}
	
	public void assembleHeatConvectSurface(double[] constConvect, BVPType tp, SRMatrix A, int... bdLabel) {
		assembleBoundaryMatrix(constConvect, FES3H81::elementConvectSurface, tp, A, bdLabel);
	}
	
	public void assembleHeatConvectSurface(double[][] vectorConvect, BVPType tp, SRMatrix A, int... bdLabel) {
		assembleBoundaryMatrix(vectorConvect, FES3H81::elementConvectSurface, tp, A, bdLabel);
	}
	
	public static void elementConvectSurface(int[] ele, double[][] xyz, double[][] h, BVPType tp, double[][] eleMatrix) {
	    double[] g = h[0];
	    double[] x = MVO.col(xyz, 0);
	    double[] y = MVO.col(xyz, 1);
	    double[] z = MVO.col(xyz, 2);
	    double xi = 0, et = 0, detJ  = 0.0, wxarea, c;
	    double[] u = new double[4];
	    MVO.fill(eleMatrix, 0.0);
	    if (tp == BVPType.CONSIST_MASS) {
	    	for(int i = 0; i < 2; i++) {
	    		xi = GP2[i];
	    		for(int j = 0; j < 2; j++) {
	    			et = GP2[j];
	    			detJ = CG.ISOMap3DTo2D(CG::Q4ShapeFunction, x, y, z, xi, et, u);
	    			c = MVO.dot_product(g, u);
	    			wxarea = GW2[i] * GW2[j] * c * detJ;
	    			for(int ii = 0; ii < 4; ii++) {
	    				for(int jj = 0; jj < 4; jj++) {
	    					eleMatrix[ii][jj] += u[ii] * u[jj] * wxarea;
	    				}
	    			}
	    		}
	    	}
	    } else if (tp == BVPType.LUMP_MASS) {
	    	for(int i = 0; i < 2; i++) {
	    		xi = GP2[i];
	    		for(int j = 0; j < 2; j++) {
	    			et = GP2[j];
	    			detJ = CG.ISOMap3DTo2D(CG::Q4ShapeFunction, x, y, z, xi, et, u);
	    			c = MVO.dot_product(g, u);
	    			wxarea = GW2[i] * GW2[j] * c * detJ;
	    			for (int si = 0; si < 4; si++) {
						eleMatrix[si][si] += u[si] * wxarea;
					}
	    		}
	    	}
	    } else
	    	throw new IllegalArgumentException("Wrong BVPType!");
	}
	

	public void computeGradient(double[] u, double[][] gradU) {
		double xi, et, zt;
		double[] uu = new double[8];
		double[] v = new double[8];
		double[] dx = new double[8];
		double[] dy = new double[8];
		double[] dz = new double[8];
		double[] count = new double[u.length];
		for (int t = 0, n = mesh.getNt(); t < n; t++) {
			int[] ele = mesh.getElement(t);
			mesh.getCoordPerElement(ele, coordRegion);
			for (int i = 0; i < 8; i++) {
				uu[i] = u[ele[i]];
			}
			double[] dxUU = new double[8];
			double[] dyUU = new double[8];
			double[] dzUU = new double[8];
			double detJ = 0.0;
			for (int i = 0; i < 2; i++) {
				xi = GP2[i];
				for (int j = 0; j < 2; j++) {
					et = GP2[j];
					for (int k = 0; k < 2; k++) {
						zt = GP2[k];
						detJ = CG.ISOMap3D(CG::H8ShapeFunction, MVO.col(coordRegion, 0), MVO.col(coordRegion, 1),
								MVO.col(coordRegion, 2), xi, et, zt, v, dx, dy, dz);
						int idx = k + i * 4 + j * 2;
						dxUU[idx] = MVO.dot_product(uu, dx) * detJ;
						dyUU[idx] = MVO.dot_product(uu, dy) * detJ;
						dzUU[idx] = MVO.dot_product(uu, dz) * detJ;
					}

				}
			}
			double[] temp = MVO.matmul(Constants.GaussMatrix3D, dxUU);
			double[] temp2 = MVO.matmul(Constants.GaussMatrix3D, dyUU);
			double[] temp3 = MVO.matmul(Constants.GaussMatrix3D, dzUU);
			for (int i = 0; i < 8; i++) {
				gradU[0][ele[i]] += temp[i];
				gradU[1][ele[i]] += temp2[i];
				gradU[2][ele[i]] += temp3[i];
				count[ele[i]] += detJ;
			}
		}

		for (int i = 0; i < mesh.getNv(); i++) {
			gradU[0][i] /= count[i];
			gradU[1][i] /= count[i];
			gradU[2][i] /= count[i];
		}
	}

	public void computeGradient2(double[] u, double[][] gradU) {
		int len = u.length;
		ParamCheck.checkEqual(len, this.ndof);
		ParamCheck.checkEqual(len, gradU[0].length);
		ParamCheck.checkEqual(len, gradU[1].length);

		double[][] coord = new double[8][3];
		double[][] grad = new double[8][3];
		double[] volumeCount = new double[this.GetNdof()];
		double[] eleU = new double[8];

		for (int i = 0, n = mesh.getNt(); i < n; i++) {
			int[] ele = mesh.getElement(i);
			mesh.getCoordPerElement(ele, coord);
			for (int j = 0; j < 8; j++) {
				eleU[j] = u[ele[j]];
			}
			CG.gradientValue3DH8(coord, eleU, grad);
			double volume = CG.volume3DH8(coord);
			for (int j = 0; j < 8; j++) {
				gradU[0][ele[j]] += grad[j][0] * volume;
				gradU[1][ele[j]] += grad[j][1] * volume;
				gradU[2][ele[j]] += grad[j][2] * volume;
				volumeCount[ele[j]] += volume;
			}
		}
		for (int i = 0, n = mesh.getNv(); i < n; i++) {
			gradU[0][i] /= volumeCount[i];
			gradU[1][i] /= volumeCount[i];
			gradU[2][i] /= volumeCount[i];
		}
	}

	public double computeL2Norm(double[] u) {
		double result = 0.0;
		double[] v = new double[8];
		double[] uu = new double[8];
		double[] dx = new double[8];
		double[] dy = new double[8];
		double[] dz = new double[8];
		double xi, et, zt, detJ, wxarea, u0;
		for (int t = 0; t < mesh.getNt(); t++) {
			int[] ele = mesh.getElement(t);
			mesh.getCoordPerElement(ele, coordRegion);
			MVO.getSlice(u, ele, uu);
			for (int i = 0; i < 2; i++) {
				xi = GP2[i];
				for (int j = 0; j < 2; j++) {
					et = GP2[j];
					for (int k = 0; k < 2; k++) {
						zt = GP2[k];
						detJ = CG.ISOMap3D(CG::H8ShapeFunction, MVO.col(coordRegion, 0), MVO.col(coordRegion, 1),
								MVO.col(coordRegion, 2), xi, et, zt, v, dx, dy, dz);
						u0 = MVO.dot_product(uu, v);
						wxarea = GW2[i] * GW2[j] * GW2[k] * detJ;
						result += u0 * u0 * wxarea;
					}

				}
			}
		}
		return Math.sqrt(result);
	}

	public double computeH1SemiNorm(double[] u) {
		double result = 0.0;
		double[] v = new double[8];
		double[] uu = new double[8];
		double[] dx = new double[8];
		double[] dy = new double[8];
		double[] dz = new double[8];
		double xi, et, zt, detJ, wxarea, dxU, dyU, dzU;
		for (int t = 0; t < mesh.getNt(); t++) {
			int[] ele = mesh.getElement(t);
			mesh.getCoordPerElement(ele, coordRegion);
			MVO.getSlice(u, ele, uu);
			for (int i = 0; i < 2; i++) {
				xi = GP2[i];
				for (int j = 0; j < 2; j++) {
					et = GP2[j];
					for (int k = 0; k < 2; k++) {
						zt = GP2[k];
						detJ = CG.ISOMap3D(CG::H8ShapeFunction, MVO.col(coordRegion, 0), MVO.col(coordRegion, 1),
								MVO.col(coordRegion, 2), xi, et, zt, v, dx, dy, dz);
						dxU = MVO.dot_product(uu, dx);
						dyU = MVO.dot_product(uu, dy);
						dzU = MVO.dot_product(uu, dz);
						wxarea = GW2[i] * GW2[j] * GW2[k] * detJ;
						result += (dxU * dxU + dyU * dyU + dzU * dzU) * wxarea;
					}

				}
			}
		}
		return Math.sqrt(result);
	}
}
